arXiv: 1509.06763v2 [quant-ph] 4Jul2016 


Practical, Reliable Error Bars in Quantum Tomography 


Philippe Faist* and Renato Renner 
Institute for Theoretical Physics, ETH Zurich, 8093 Switzerland 
(Dated: July 5, 2016) 

Precise characterization of quantum devices is usually achieved with quantum tomography. However, most 
methods which are currently widely used in experiments, such as maximum likelihood estimation, lack a well- 
justified error analysis. Promising recent methods based on confidence regions are difficult to apply in practice 
or yield error bars which are unnecessarily large. Here, we propose a practical yet robust method for obtaining 
error bars. We do so by introducing a novel representation of the output of the tomography procedure, the 
quantum error bars. This representation is (i) concise, being given in terms of few parameters, (ii) intuitive, 
providing a fair idea of the “spread” of the error, and (iii) useful, containing the necessary information for 
constructing confidence regions. The statements resulting from our method are formulated in terms of a figure 
of merit, such as the fidelity to a reference state. We present an algorithm for computing this representation 
and provide ready-to-use software. Our procedure is applied to actual experimental data obtained from two 
superconducting qubits in an entangled state, demonstrating the applicability of our method. 


Introduction .—Recent experimental developments have 
demonstrated increasingly precise manipulation and control 
of quantum systems, paving the way towards the hopeful im¬ 
plementation of a quantum computer [1-12]. The successful 
outcome of an experiment is usually certified using quantum 
tomography. This is the task of inferring the quantum state 
of a device from statistics of measurements on many copies of 
the system [13-19]. Several methods perform this task and are 
widely used, such as maximum likelihood estimation [20, 21]. 

In the realistic regime where finite data are collected, the er¬ 
ror bars provided by most methods which are widely applied 
in current experiments [19, 22-24] are typically ill justified 
and may lead to deceiving conclusions [25-27]. To remedy 
this problem, Blume-Kohout [27] and Christandl and Ren¬ 
ner [28] resort to confidence regions. These are regions in 
state space of all density matrices in which the state lies with 
high probability. In contrast to Bayesian methods [25], the 
reliability statements do not depend on any prior distribution. 
However, confidence regions are a priori difficult to construct 
explicitly [29]. Furthermore, they are designed for worst-case 
scenarios and are often not representative of the intuitive ex¬ 
tent of the error. 

Our main result is a novel representation of the output of 
the tomography procedure—a summary of what the tomo¬ 
graphic data tells us about the state of the system—which we 
call quantum error bars. This description is (i) concise, being 
given in terms of a few parameters only, (ii) intuitive, provid¬ 
ing a fair idea of the “spread” of the error, and (iii) useful for 
precise statements, containing all necessary information for 
constructing confidence regions. Our method, in particular, 
inherits the mathematical robustness of the confidence region 
approach. 

The quantum error bars are designed to mimic the role of 
classical error bars. Classically, an error bar typically repre¬ 
sents the standard deviation of the distribution of a physical 
quantity caused by noise or statistical errors; this distribution 
is usually assumed to be Gaussian. Observe that, precisely, 
classical error bars (i) are a concise description of the error, 
(ii) provide a fair, intuitive idea of the “spread” of the quantity 
of interest, and (iii) allow us to calculate precise statements 
such as the required error interval to consider (e.g., 5 standard 


deviations) for a specific requested certainty level (e.g., one in 
a million). 

Our statements are formulated in terms of a figure of merit 
which can be chosen freely. Our method works best when 
the figure of merit is the fidelity to a pure target state, the 
expectation value of an observable, or the trace distance to any 
reference state. This encompasses most tomography settings. 

The quantum error bars are constructed as follows. The 
input is the experimental data from a general quantum tomog¬ 
raphy experiment. Then we construct a particular distribution 
jj. (/) of the chosen figure of merit /, which has the property of 
containing the necessary information to construct confidence 
regions at any confidence level using the method of Ref. [28]. 
We show that in a wide range of situations and for a class of 
figures of merit, the distribution /j(/) can be approximated 
by a simple analytical expression with three parameters. The 
quantum error bars are then straightforwardly deduced from 
these parameters. 

We provide a simple numerical algorithm to obtain the 
quantum error bars from the measurement data. By fitting 
a numerical approximation to /j(/) with our approximate an¬ 
alytical model, we obtain the values of the parameters of the 
model which directly translate to the quantum error bars. The 
practicality of our method is demonstrated by applying it to 
experimental data from two superconducting qubits. 

Our work complements a vast literature which has pro¬ 
vided error analyses for experiments [30—42] as well as ex¬ 
plicit schemes [43-52], by introducing the novel concept of 
quantum error bars. The complexity of such schemes have 
also been investigated [53, 54] and numerical techniques put 
forward [25, 55-57]. Furthermore, a number of contributions 
propose measurement schemes for fidelity estimation [58, 59], 
tomography of matrix product states [60], estimation of low- 
rank states [61, 62], and permutationally invariant tomogra¬ 
phy [63-65]. An experiment following such schemes would 
achieve target benchmarks more efficiently, and it could still 
be analyzed using our procedure, the latter being applicable to 
any measurements. 

The rest of this Letter is structured as follows. First, we 
briefly explain our quantum tomography setup and the con¬ 
cept of confidence regions. We then derive our main technical 
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FIG. 1. Setup of quantum tomography. Measurements are taken on n 
copies of a quantum system. The outcomes allow us to infer what the 
state of the quantum system is. In this example a qubit is measured 
using Pauli operators. Here, the experimental data are most consis¬ 
tent with the state being |t), located at the top of the Bloch sphere 
(in green). However, because only finite data are collected, there 
is an uncertainty associated with this statement. In the method of 
Ref. [28], a distribution (p) (the red gradient) is determined from 
the data, from which confidence regions can be constructed (delim¬ 
ited by the dotted line). These are regions in state space in which the 
state lies with high probability. 


results, namely, the definition of its approximate theo¬ 
retical model, and the algorithm to estimate fr(/) numerically. 
Finally, we demonstrate the applicability of our method on ex¬ 
perimental data. 

Quantum Tomography Setup .—A large number n of copies 
of a quantum system are measured using independent, possi¬ 
bly different, measurement settings (Fig. 1) [66]. We list all 
the of the distinct positive operator-valued measure (POVM) 
effects in one set {E;^}, and denote by ni^ the number of times 
the POVM effect Ej^ was observed. We then construct the like¬ 
lihood function, which will be needed in our analysis. It is 
defined as the probability with which the observed data would 
occur if the true state were « copies of p, 

A(p) = Pr[observed data | p] = ]~[(tr[£'(.p])”* , (1) 

k 

along with the log-likelihood, 

-^(P) =-21nA(p) =-2^nA.lntr(£'<.p) , (2) 

k 

with a conventional (—2) factor [27, 33]. 

Confidence Regions .—In the following, we briefly review 
the method of Ref. [28] for constructing confidence regions, 
on which our method is based. 

Confidence regions of confidence level 1 — e are defined 
as regions in state space which contain the true state with 
a probability of at least 1 — e. Crucially, it is the complete 
procedure of assigning a region to tomographic data which 
is certified and not the particular region itself (despite the 
slightly misleading terminology). More precisely, for a par¬ 
ticular “true” state Ptrue^ the measurement outcomes observed 
in the tomography procedure are only one possible outcome 
data set among the enormous amount of theoretically possi¬ 
ble data sets. Now, a data analysis procedure associates with 
each observed data set a corresponding region in state space. 
This tomography procedure is said to yield confidence regions 


of confidence level 1 — e if, for any true state Ptme, the to¬ 
mography procedure associates with the observed data set a 
region which contains Ptme. except for some data sets with to¬ 
tal probability £. In other words, the complete tomography 
procedure is successful except with probability e, in which 
case the observed data set may cause the procedure to report a 
bad region. These “exceptional data sets” may be interpreted 
as misleading but highly unlikely situations. For example, if 
we flip a fair coin many times and observe the sequence of 
all “heads,” any reasonable inference scheme would wrongly 
report that the coin is highly biased. However this outcome 
only happens with disproportionately small probability; intro¬ 
ducing the parameter e above allows us to disregard such ex¬ 
tremely unlikely cases. 

The method of Ref. [28] is formulated using the estimate 
density Pb" [67], defined as 

Ab«(p) = — A(P) , (3) 

Cb” 

where cb” is a normalizing factor such that / dp Pb"{p) = T 
and where dp is the Hilbert-Schmidt measure normalized 
such that f dp = 1 [68, 69]. The main result of Ref. [28] is 
a criterion for certifying a procedure for yielding confidence 
regions of confidence level 1 — e. The criterion is the fol¬ 
lowing: the procedure should map to any tomographic data 
(essentially) a region R in state space which satisfies 

PB"iP)dp = l -, ( 4 ) 

Jr poly(n) 

i.e., which has high weight under the distribution ps" [70]. 

Confidence Regions for a Figure of Merit. —We may now 
use this criterion to devise an explicit procedure for construct¬ 
ing confidence regions, where the regions R are chosen to be 
defined via level sets of a figure of merit. 

A figure of merit /(p) may be any function of the quantum 
state. For example, /(p) = F^{p, | Vfiief)(Vfitef|) expresses the 
fidelity to a reference state | Vfitef) ■ The reduced distribution of 
the estimate density Pb" (p) onto the figure of merit / is given 
by 

P(/) = JdpPB"{p)5{f{p)-f) , (5) 

where 5() denotes the Dirac delta function. 

Now, fix a threshold value /, and consider the region Rf 
in state space consisting of all states whose figure of merit is 
greater than or equal to / (Fig. 2). The weight of the region 
Rf according to the distribution Pb"{p) is exactly given by 
p(/') df'. Inverting this reasoning, for any e, we can 
find the maximum threshold value / required for a region R f 
to encompass a particular weight 1 — e/poly (n); we know that 
this region is essentially a confidence region by the criterion 
of Ref. [28]. (If the figure of merit is such that smaller values 
of /(p) are desirable, such as the trace distance to a reference 
state, then Rf is defined with / as an upper, rather than lower, 
threshold value.) 

We arrive at a first important observation: if we find a sim¬ 
ple characterization of the function p{f), then we are capa¬ 
ble of constructing confidence regions in terms of / for any 
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FIG. 2. Construction of confidence regions from the distribution 
/I (/) on the figure of merit. High weight intervals with respect to 
/t (/) (left plot) correspond to high weight regions in state space with 
respect to flgK (p) (right diagram) which are (essentially) confidence 
regions, according to Ref. [28]. 


confidence level (See Appendix D for how to transpose the 
5-enlargetnent in [28] into a shift of the threshold value /.). 

Determining ll{f) Numerically . — We propose a practical 
procedure which determines a numerical estimate of /f(/). 
We resort to a Monte Carlo-type scheme known as the 
Metropolis-Hastings algorithm [71] (cf. also Refs. [72, 73]). 
This algorithm is a standard, well-tested scheme widely used 
in computational physics—for instance, to simulate the be¬ 
havior of statistical systems at finite temperature [74] —and 
there are standard methods for controlling the uncertainties 
resulting from the use of this procedure [75]. Using this algo¬ 
rithm, we conduct a random walk in the quantum state space 
and produce random samples distributed according to the dis¬ 
tribution jLIb'i (p). By collecting the values of /(p) at the sam¬ 
pled points into a histogram, we obtain an estimate for p(/). 
(See Appendix A for the details of the random walk proce¬ 
dure). 

Theoretical Model for —It turns out that, for a se¬ 

lection of common figures of merit, we may understand the 
numerical estimate of p(/) with a theoretical model. Sup¬ 
pose /(p) is the fidelity to a pure reference state, the expecta¬ 
tion value of an observable, or the trace distance to any refer¬ 
ence state. Then, under some reasonable assumptions [76], we 
derive the following approximate theoretical model for p (/) 
(see Appendix B), 

M (/)« C (/-/!)'" • (/-'") , (6) 

with three fit parameters ai, a 2 , and m (with m ^ 0 ), and one 
constant normalization factor C', h is a constant depending 
only on the choice of the figure of merit. Specific values of 
the constant li for some figures of merit are summarized in 
Table I. 

The parameters (a 2 ,ai,m) are then mapped onto new pa¬ 
rameters which are more representative of the shape of the 
function. The latter is viewed as a “skewed” Gaussian (see 
Appendix C). The parameter /o determines the position of 
the peak, the parameter A is the half width of the “deskewed” 
Gaussian, and 7 characterizes the deviation from a perfect 
Gaussian. The parameters (/o,A, 7 ) are the quantum error 
bars. 

Application to Experimental Data . — We have applied the 


Inp (/) ~ —a2X^ — aix + m\ 
Figure of merit / (p) 

nx + c, where: 

X = 

(p, 1 V4;ef)( V^Ref 1) = ( V^Ref 1 P 1 V^Ref) 
U(p,PRef) = slip-PRef 111 

Observable (A)p 

1 -/ 

/ 

a-foif-a 


TABLE I. Theoretical fit model for some selected figures of merit. 
Here | Vfiief) denotes any pure state, and p^gf any pure or mixed state. 
We use the notation D (p, cr) for the trace distance and (A) p = tr (Ap) 
for the expectation value of an observable A. The value a is an ex¬ 
tremal value of (A) p for valid density matrices p close to the region 
of interest, and x should be chosen as x = a — f (resp. x = f — a) if a 
is a maximal value (resp. minimal value). If the extremum point of A 
is far from the region of interest, the logarithm term in the model can 
be dropped, as the exponential will dominate the volume term, and a 
can be absorbed into the other factors. 
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FIG. 3. Analysis of measurement data from two superconducting 
qubits prepared in a Bell state. We determined effective measure¬ 
ment operators which model the noise in the measurement process. 
The histogram of the fidelity to the target state |) 7 ) (the blue data 
points), produced using our procedure, fits well to our theoretical 
model in Table I. The quantum error bars are a concise, intuitive, 
and precise characterization of the fit model, which is interpreted as 
a skewed Gaussian function. The parameter /q is the peak maximum, 
A is the half width of the original Gaussian, and 7 characterizes the 
skewing in terms of the displacement of the sides of the peak from 
the exact Gaussian at the relative height 1 /e. This example involving 
experimental data demonstrates a good level of practical applicability 
of our method. 


algorithm to experimental data from two superconducting 
qubits prepared in a Bell state according to the setup de¬ 
scribed in Refs. [10, 77]. The data were kindly provided by 
the authors of Ref. [10]. The two qubits were measured us¬ 
ing slightly noisy individual Pauli operators, with a total of 
n — 55677 measurements. The numerical estimation of p. (/) 
corresponding to the fidelity to the target Bell state is depicted 
in Fig. 3. (See Appendix F for details of the analysis of the 
experiment, including the modeling of the measurements [78] 
into effective POVM operators.) 

Quantum Error Bars .—The quantum error bars (/o. A, 7 ) 
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displayed in Fig. 3 are a concise and useful description of 
the error analysis, from which reliable operational statements 
can be made. Indeed, they provide the necessary information 
for constructing confidence regions for any given confidence 
level. 

Also, as seen in Fig. 3, our error bars have the intuitive in¬ 
terpretation as representing the “spread” of the figure of merit 
according to /r (/). As such, the error bars are much smaller 
than the size of a confidence region for a small epsilon in a 
worst-case scenario, and they are in fact of comparable size 
to those obtained by bootstrapping [22, 24, 27, 41, 79] (see 
Appendix G). 

Discussion .—Our work bridges the apparent gap between 
carrying out a mathematically rigorous, well-justified error 
analysis and using an ad hoc procedure yielding smaller error 
bars. The quantum error bars provide a convenient and precise 
representation of the information provided by the tomography 
procedure. 

While the fit model for /r(/) is subject to some assump¬ 
tions and approximations, it applies well to many examples 
studied by the authors in developing this work—for n ^ 100 
total measurements already—and has been tested with up to 
five qubits. Note that the numerical procedure is not subject 
to these assumptions, and a deviation from the fit model could 
easily be noticed in some extreme examples considered (for 
example, with goodness-of-fit measures). A further detailed 
discussion on the reliability of our method is presented in Ap¬ 
pendix H. 

It is relatively straightforward to apply our method to exper¬ 
imental setups consisting of a few qubits. Our procedure is re¬ 
stricted neither to particular measurement settings nor to spe¬ 
cific quantum states, and it applies, for example, to adaptive 
tomography. In general, noise in the measurement procedure 
has to be modeled into effective POVM effects analogously 
to our approach for the two superconducting qubits. (In con¬ 
trast, other approaches do not require this [80-82].) We have 
developed a software which implements our procedure [83] 
that is expected to be directly applicable to most experimental 
settings. 

For worst-case scenarios such as quantum cryptogra¬ 
phy [84], it is still desirable to improve the methods for ex¬ 
plicitly constructing confidence regions. We do anticipate that 
the bounds used in Ref. [28] may be tightened to yield smaller 
confidence regions for the same confidence level. If the con¬ 
struction is not altered, the procedure presented here would 
not require any change, as the same histograms may still serve 
for constructing confidence regions using the tightened proof. 

We also insist that our results do not rely on any particular 
interpretation of “probability,” such as a Bayesian or a fre- 
quentist one. This is because we consider experiments which 
can, in principle, be repeated arbitrarily many times, which 
is a regime where these interpretations are equivalent [28]. 
Nonetheless, the Bayesian viewpoint is convenient, as the dis¬ 
tribution jj. (/) happens to coincide with the Bayesian poste¬ 
rior corresponding to an agent starting the tomography proce¬ 
dure with a Hilbert-Schmidt uniform prior. 

Furthermore, even though our results are formulated in the 
context of quantum state tomography, the same procedure 


may be applied to quantum process tomography [85, 86]. In¬ 
deed, the Choi-Jamiolkowski isomorphism implies that deter¬ 
mining a quantum process is mathematically the same as de¬ 
termining a bipartite quantum state. 
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SUPPLEMENTAL MATERIAL 


In this appendix, we provide a detailed description of how 
our method is implemented, how it is applied to practical ex¬ 
amples, as well as additional discussions referred to from the 
main text. An overview of our method is given in Fig. 4. 

Our software, along with instructions for download and 
use, may be downloaded at the location: https://tomographer. 
github.io/tomographer. 


Appendix A: Procedure for determining /x(/) using a 
Metropolis-Hastings random walk 

The figure of merit is given as a function / (p) of the quan¬ 
tum state. For example, the squared fidelity to a reference 
state pRef is represented as f[p)=F^ (p,pRef). 

Recall that the relevant object in the method of Christandl 
and Renner is the estimate density (p), given by Eq. (3) of 
the main text. 

Given the figure of merit /(p) of interest, the reduced dis¬ 
tribution on this of merit of ps" (p) is 

M(/) = JdpPB'’{p)d[f{p)-f] , (Al) 

where the Dirac delta ensures the integration is performed 
over the shell of states in state space which have the given 
figure of merit /. The quantity p (/) corresponds to the total 
weight given by ps” (p) to all states with a given fixed figure 
of merit /. 

In the following, we develop a method to compute p (/) 
numerically. We resort to a Monte Carlo-type scheme known 
as the Metropolis-Hastings algorithm [71] (cf. also Refs. [72, 
73]). This scheme is widely used in computational physics, 
for instance to simulate the behavior of statistical systems at 
finite temperature [74]. This algorithm conducts a random 
walk which produces random samples distributed according to 
a given distribution P (x). The parameters of the algorithm are 
a starting point xq as as well as a “jump distribution” Q (x'|x). 
The jump distribution is assumed to be symmetric {Q (x'|x) = 
Q (x|x')), and is used to update the current step in the random 
walk. (For example, Q {x' |x) is often chosen as a Gaussian in 
some relevant coordinates centered at x). The i-th step of the 
random walk goes as follows: 

1. Choose a new candidate point x' according to Q (x'|x,); 

2. Calculate a = P{x!)/P{xi). If a > 1, then set x,+i := x' 
unconditionally; if a < 1 , then decide randomly to set 
Xj+i := x' with probability a, or else to set x,+i := x;. 

The sequence of points {x, }, albeit correlated, are then asymp¬ 
totically distributed according to the distribution P(x). 




Measurement Data 


> tomo 


Numerical procedure: 
Metropolis-Hastings Random Walk 






Error Bars 


FIG. 4. Overview of our quantum tomography analysis and how 
to apply our method. The measurement data are the input to our 
procedure. Our numerical method, which is based on a Metropolis- 
Hastings random walk in state space, outputs a histogram of a chosen 
figure of merit. We provide software accomplishing this [83]. This 
histogram is a numerical approximation to the distribution p (/) of 
the figure of merit. In a second step, this numerical estimate is fitted 
by a theoretical model. The fit can be done, for instance, using MAT- 
LAB. This yields a full description of the relevant distribution which 
allows to construct in principle confidence regions for any confidence 
level. This description is given in terms of three parameters, which 
are then effectively the “error bars.” 


In order to calculate the quantity p (/), we draw a large 
number of random samples in the quantum state space accord¬ 
ing to the distribution ps'’ (p), i.e. with ps" (p) playing the 
role of P (x). The histogram of values / (p) evaluated at those 
samples then provide a numerical estimate of p (/). Crucially, 
it is not necessary to calculate the normalization constant cb" 
in the definition of Pb” (Eq. (3) of the main text) because in 
the Metropolis-Hastings algorithm we only have to evaluate 
ratios of probabilities. 

For our random walk, we represent a quantum state p of 
dimension li by a square complex matrix T with trTT^ = 1, 
such that p = rr^. To any such T corresponds a valid density 
matrix p, and to any density matrix p corresponds at least one 
such T (e.g. T — p '/^). Additionally, the constraint txTT'^ = 1 
corresponds to requiring that the components of T, real and 
imaginary parts taken separately into a real vector y, lie on the 
surface of the unit {2d^ — l)-hypersphere. Random density 
matrices may be sampled from the Hilbert-Schmidt measure 
by choosing such random points on this hypersphere [68]. In 
fact, the matrix entries 7]y of T are simply the components 
of a vector jy/) of dimension d^ which purifies p. Indeed, 
if we trace out the second system from lyr) = LiyT/flOA ® 
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1;%, we obtain ItbIyXyI = Liji' TijT;^j\i){i'\ = = p. We 

hence choose to perform a Metropolis-Hastings random walk 
on the {2d^ — 1 )-hypersphere corresponding to the possible 
T matrices. The jump candidate is calculated from a point 
yi by choosing a vector cu of Icfi normally distributed values 
and setting f = (y; -f Tjstep®)/||yi + I7,step®||, where T7step is a 
chosen step size. The jump distribution obtained in this way 
is symmetric. 

We follow the the prescriptions given in Ref. [75] for the 
correct usage and appropriate error analysis of the Metropolis- 
Hastings algorithm. First, since the random starting point is 
likely to be a point which has very little weight under P{x), 
the random walk needs to equilibrate, or thermalize, until it 
reaches points which have a non-negligible values of P{x). 
This first set of points traversed until the walk has thermalized 
should be discarded. Also, because consecutively collected 
samples may be very correlated, it is useful to keep only one 
sample in a certain number A^sweep (the “sweep size”), while 
throwing away each time the A^sweep — 1 points between two 
recorded samples. In our examples, the sweep size A^sweep is 
chosen of the order of 1 / rjstep', this gives at least the chance 
to the random walk to traverse all of state space between two 
recorded samples. Errors on the final numerical histogram 
points may be estimated either by calculating the standard de¬ 
viation of independent runs of the simulation, or by a binning 
analysis which takes into account the correlation of the sam¬ 
ples during a single run. We refer to Ref. [75] for a detailed 
discussion of these techniques. 

In our case, for each histogram bin, we associate to each 
recorded sample the value 1 if the point is in the bin, or 0 
otherwise. The final numerical estimate of p (/) is produced 
by averaging the time series for each bin, which corresponds 
up to a constant to generating a histogram of counts. These 
time series are suitable for use in a binning analysis to obtain 
error bars on the numerical estimate of p (/). 

Appendix B: The fit model for p (/) 

We now derive an approximate theoretical model to fit our 
numerical histogram. This is useful in order to provide a suc- 
cint description of the result in terms of only a few parame¬ 
ters. It also serves as a consistency check allowing us to assert 
that our results are well understood from a theoretical point of 
view. 

In general, the function p (/) can be very complicated, so 
an exact analytical description is unlikely. Rather, our goal 
is to find a decent approximation of this function in a region 
close to where ps" has high weight. 

In fact, the bell shape of the curves in Fig. 3 of the main text 
is typical when the figure of merit is taken to be the fidelity to 
a pure target state, the expectation value of an observable or 
the trace distance to any reference state. Intuitively, this shape 
is the result of the concurrence of two effects: a volume factor 
reflecting the increasing surface of a shell of fixed figure of 
merit as we get far from the reference point, and the approx¬ 
imately exponential decrease of the likelihood function itself. 
For example, in the case of the trace distance to the maximum 


likelihood estimate, there are many more states with high dis¬ 
tance to jOMLE than there are very close—this is the increasing 
volume factor. The following derivation makes this argument 
more precise. 

Let’s now derive the fit model. We parameterize p with 
a generalized Bloch vector [87, 88 ]. Take an orthonor¬ 
mal basis {Aj} of the Lie algebra su{d), with j = 
and M = d^ — 1. The Aj are hermitian, traceless and obey 
trAjAji = 5jji (an example are the normalized generalized 
Gell-Mann matrices [87, 88 ], or, for k qubits, the normalized 
tensor product of Pauli operators). We may now write a gen¬ 
eral state p as p {a j) = {l/d)t+Y,jOjAj with real coefficients 
aj obeying some nontrivial constraints such that p > 0. The 
Hilbert-Schmidt distance is given by the Euclidean distance 
of the generalized Bloch vectors, tr[(^p{aj) — p{bj))^'\ = 

Now because the Hilbert-Schmidt measure dp 
is induced by the Hilbert-Schmidt metric [ 68 , 69], we may 
write (Al) as 

M (/) = ^ /5 [f{aj)-f] , (Bl) 

with a new constant c' and A (p) = — 21nA (p), and where im¬ 
plicitly the arguments to A (•) and / (•) are to be transformed 
into p appropriately. 

The conditions that the aj have to satisfy in order to rep¬ 
resent a positive semidefinite matrix are complicated [87, 89, 
90]. However, it turns out that we don’t need to know the 
exact form of these constraints. Rather, we assume that: 

(i) / has an extremal value close to the region of interest 
(viz., near Pmle); 

(ii) the surface of a shell of states of a given figure of merit 
/ tends to zero as / tends to this extremal value. 

These assumptions are rather natural and are indeed automat¬ 
ically satisfied if / (p) is one of the cases considered in the 
main text (the squared fidelity to a pure reference state, the 
trace distance to a reference state PRef, or the expectation 
value of an observable). In the case of a distance measure, 
such as the trace distance, the extremum is usually zero at the 
reference state itself, and the surface of the shell of states with 
very small distance to pRef clearly shrinks to zero. In the case 
of the expectation value of an observable, the extremum is at¬ 
tained at the border of state space. Because the border of state 
space has no flat facets, the surface of a shell of given expec¬ 
tation value (a hyperplane intersected with state space) also 
shrinks to zero as we approach the border. Furthermore recall 
that the squared fidelity to a pure reference state | y/Ref) can be 
written as the expectation value of the observable | V 6 tef)(V^Ref|- 

Denote by pRef a relevant reference point where the figure 
of merit is extremal, and let such that pRef = (1 /d)l -f 
(We recycle the notation pRef since whenever the 
figure of merit is a distance measure to a reference state, the 
same reference state is to be used here.) We go to hyperspher- 
ical coordinates (r,n) centered at the reference point 
with d^Qj = drdQ.r^^\ and introduce the change of vari- 
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ables r ^ r' = f{r,D.): 
1 


M (/) = ^ /5 [/(r,n) -/] 
dr'da 


1 


pw-l 

df 

-\ 


dr 



da 


yM-\ 

df 

-f 


dr 



5[,J -f] 
e-2^(/.n) ^ (32) 


where in the last two integrals the terms in square brackets 
are to be evaluated at the points r which satisfy r' = /(r, a) 
and / = f{r,a), respectively. Note that the figure of merit 
f{r,a) must be invertible for fixed a and for r ^ 0; this is 
usually the case with our choice of above. Note also that 
Expression (B2) is in fact still exact, albeit very difficult to 
calculate explicitly. To proceed further, we will use Laplace’s 
approximation, and assume that the main contribution to the 
integral is a region close to a single point Ho on the shell of 
fixed figure of merit / where the integrand is maximal. Then, 
we have 


(B2) 


1 

M-l 

df 

-f 


r 

dr 





g-^A(/,ao) ^ 


(B3) 


where w [f, Ho] is a “width factor,” which accounts for the 
total weight of the peak in Laplace’s approximation, including 
a possible truncation of the peak caused by the border of state 
space. Note that the condition (ii) above for the figure of merit 
ensures that w [/jfio] doesn’t explode when / approaches the 
extremum of / (p). 

At this point, we need to make further assumptions about 
the behavior in / of the individual terms in (B3). Lirst, we 
consider the regime in which the likelihood is close to Gaus¬ 
sian in the Hilbert-Schmidt coordinates. 


A (5) « Ao + Aa • (fl - + id- d^^^fXe {a - , 

(B4) 


(we shift fl by without loss of generality and for later prac¬ 

tical reasons; also Ag is a symmetric matrix). This is generi- 
cally the case for most practical scenarios where a reasonable 
amount of measurements were taken. 

Second, we need to assume something about the figure of 
merit /(r,n): we’ll suppose that 

fir,a) = rgia)+h, (B5) 


where h is some known constant, and g{a) some function. 
The figures of merit considered in the main text automat¬ 
ically satisfy this assumption. Lirst, if the figure of merit 
is any distance measure to pRef which is given by a norm, 
such as the trace distance, then /(r,fi) = ||p (r,fi) — PRef|| = 
\\Y.j{aj{r,a) - a^®f) Aj\\ = r||£yf2;Ay||, recalling that our hy- 

perspherical coordinates are defined by d{r,a) = a^®^-|-rn 
with a the unit vector in the direction a. Also, /(p) obeys 
property (B5) if it is the expectation value of an observable, 
/(p) = tr(pVT); with a = we can write /(p) = 

tr([(l/£/) it) =rn-vv-|-(a^®^-w-|-trlT/t/), where 
w is the vector with components wj = tr(AylT). Recall in 
this case that is on the border of state space. Lurther- 
more, recall that the squared fidelity to a pure reference state 
can be written as the expectation value of an observable, 
lrRef)(rRef|) = tr(p|V/Ref)(rRefI)- However, if the fig¬ 
ure of merit is the fidelity or purified distance to a mixed ref¬ 
erence state, it does not satisfy in general the Ansatz (B5). 

Armed with both assumptions (B4) and (B5), we see that 
df /dr = g(a) as well as r = (/ — h) jg (fi), and we obtain 


P(/) 




■w[f,ao] 


^ 2 [T)+^A4'^0+P^o 


(B6) 


where c" = c' ■ sign(g(Ho)). The term in the exponential 
in (B6), being quadratic in r, is then also quadratic in f — h. 
At this point, we further assume that Hq (where A {f,ao) is 
minimal at fixed /) is approximately constant in /, and that 
the term w [/, fio] is either approximately constant or, being a 
volume factor, varies as a power of r, and thus of f — h. We 
finally obtain our fit model. 


P (/) « C (/ - /z)'" • (/-*) , (B7) 


with 3 fit parameters ai, 02 , m and one constant normalization 
factor C. The value m includes the (M — 1) power plus any 
contribution from the weight factor w [/, H]. The expression 
for the logarithm of p (/) is numerically more suitable for 
fitting. We thus obtain our fit model for In p (/), 


lnp(/) 


— 02 if — h)^ — a\ if — h)+m\nif — h)+c, or 

— 02 ih — f)^ — a\ ih — f) +m In (/; — /) -f c , 


(B8a) 

(B8b) 


depending on whether f ^ h for all valid / or / ^ /z for all valid /. Indeed, either of these two conditions hold as we have 
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chosen the center of our hyperspherical coordinates as an 
extremal point of /(p). Table I of the main text summarizes 
the appropriate fit model for a selection of figure of merits. 

It is further worth mentioning that for larger /, the expo¬ 
nential will dominate all the other terms; for example, in this 
regime, the details of the function w (/, H) is not relevant for 
most figures of merit. 

There are certain situations in which our approximate fit 
model fails to accurately describe the behavior of p (/). If too 
few measurements are taken, the likelihood function is not 
approximately Gaussian as we have assumed (however this is 
usually the case already for, e.g., n ^ 100 total measurements). 
Our derivation also no longer applies if Oq happens to not be 
constant with /, or if a different figure of merit is considered 
such as the fidelity to a mixed reference state. Furthermore, 
in some cases the boundary of state space might interfere with 
our approximation (it might for example constrain Ho caus¬ 
ing it to vary with /), or the Laplace method might not be 
a good approximation if A has e.g. several minima for 

fixed /. However in examples we have studied these cases al¬ 
ways caused our model to fit poorly to the numerical estimate 
in the region of the peak; we consider it very unlikely that the 
fit model fits the peak well but fails to describe the tail accu¬ 
rately. See also Appendix H for a more general discussion of 
the reliability of our method. 


Appendix C: The quantum error bars 


Now we proceed to transform the parameters {a2,ai,m) 
into more meaningful quantities, corresponding to distinctive 
features of the corresponding function. Consider the function 

y(x) =—a2X^— aix + mln(x)+c , (Cl) 

which for X > 0 exhibits a characteristic skewed bell curve as 
depicted in Fig. 3 of the main text. The function y(x) is ex¬ 
actly the model for ln/r(/), with x = s{f — h) for a constant 
h and for a sign i = ±1 depending on the figure of merit, as 
given by (B 8 ). The function (Cl) can be seen as a skewed ver¬ 
sion of a parabola with summit at (xo,y(xo)) (see Figure 5a). 
The de-skewing operation applied to y(x) consists in finding 
the parabola ydeskewed (■^) = —a{x — xq)^ +yo with the same 
the peak position and curvature as y{x). In other words, we 
find a parabola ydeskewed(^) such that the functions y{x) and 
ydeskewed(-*) must match at X = xq to second order. The maxi¬ 
mum of y(x) is at the point xq satisfying 0 = {dy/dx)\^^, that 
is, 0 = —2 02 X 0 — ai -l-m/xo. Solving for xq while ensuring 
that Xq > 0 yields 


XQ =- (—ai + \/+ Sa2m] . (C2) 

402 V ^ ) 

The condition {d^yjdx^^f;, = (t/^ydeskewed/<^.*^)|;to yields a = 
02 -ym/{ 2 x^). In terms of (o,xo,yo)5 the original parameters 




FIG. 5. Model function for p(/) as skewed Gaussian, a. The func¬ 
tion y(x) = —a 2 X^ — aix + m\nx + c (used to model lnp(/)) can be 
seen as a skewed parabola, with a “skewing operation” parametrized 
by y. The green curve Vdeskewed)-'^) is an actual parabola, b. The 
model function y(x) is fully characterized by the location of the peak 
XQ, the half width A of the “de-skewed” parabola at relative height 
1 /e, and a parameter y characterizing the shift of the sides of the 
peak at relative height 1 /e. c. The variable / relates to x via a simple 
shift and possible reflection, given by x = i- (/ — h) for a constant h 
and 5 = ±1 depending on the figure of merit (see Table 1 of the main 
text). Hence, with /o = sxq + Ii, the curve fl{f) is fully characterized 
by the parameters (/o, A, y), which we call the quantum error bars. 


{a 2 ,ai,c) read 


02 = 0 - ; 

2 x^ 

(C3a) 

2m 

fli = - 2 flXo ; 

(C3b) 

XQ 


c = yo +02X0 + 01 XQ — mlnxo . 

(C3c) 


We can already define A, which is the first of our quantum 
error bars. It is defined as 


A = 





(C4) 


The parameter A is the half width of the Gaussian function 
gWcskewedlx) relative height 1 je (Figure 5b): Indeed, the stan¬ 
dard deviation of a Gaussian is precisely the half width of the 
Gaussian peak at relative height 1/e with respect to the Gaus¬ 
sian peak maximum. In our case, A is interpreted as the half 
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width of the Gaussian, before the skewing operation is ap¬ 
plied. 


It remains to understand the effect of the m parameter in 
terms of skewing. Consider the intercepts of ydeskewed(-^) 
with the line y = yo — which are at x± = xo± 

(These points correspond to the cross-section of the peak of 
e^dcskewed(j:) at a relative height e^^.) If we view the function 
y(.r) as the result of skewing 3 'deskewed(-^) via the transforma¬ 
tion above parametrized by m, then the intercepts with the line 
y = yo — ^ are shifted by some which vary as a function of 
m. Let us determine 5x± to first order in m. For infinitesimal 
m, the equation y{x ±) = yo — ^ defining x± varies correspond¬ 
ingly as y{x± + 5x ±) -f 5y (x± -f dx±) =yo — ^. Keeping only 


the terms of first order in m we obtain 


dy 

dx 


-C± 


5x± -f 5y(;ic±) = 0 . 


(C5) 


Noting that we only need {dy/dx)\x^ to zeroth order in m, we 
have 


= — 2.(i2^± — U] = —-j- 2flXo 

x± ,m=0 

= (C 6 ) 

Also, with 5a2 = —mf{2xQ), 5ai = 2m/xQ and 5c = 
{5a2)x^ -f 5a 1 — mlnxQ, we have 


di 

dx 


5y(x±) = —( 5 a 2 )x^ — {5ai)x± -f m ln(x±) + 5c 
= ( 5 a 2 )(xo —x|) -I- (5ai)(xo —x±) -l-OT In 


x± 

xo 


= -^- (t 2 xo^ ^/^A - (<« T — • (<? +m\ni\ + -— 

Xq \ Xq j 




xo 


2 x^ 

ZXq 


mini 1 ± 




xo 


Developing the logarithm as a Taylor series in A, the first two orders cancel and we have 

±m(^'/^A)^ (—1)(=F1)^^^/^A^ 


{Cl)'. 


2x1 


4x4 


k-x^ 


Then we obtain from (C5), also recalling that a = A 


-2 


= -(¥ 2 « 


■ ^y\x± 

) y X..3 


(-1)(=F1)^^^/^A^ 


3x: 


k-x^ 


-f • 


3 T A h '' ' H 


6 x; 


0 


8x4 


2 k -Xq 


(C7) 


(C 8 ) 


(C9) 


We now introduce the skewing factor j as 


More precisely, the shift for infinitesimal m is given by 


inA^ 
^ 6x3 


(CIO) 


such that to lowest order in A, the shift of the “sides of the 
peak” given by x± for a relative height is directly propor¬ 
tional to 7 (see also Figure 5b): 


5x±«^ 7. (Cll) 


5x± = 7 - 



3 ^3/2A 
4x0 


+ ...+ 


. 3.|S{^72)+1a^' 

(^' + 3)-4' 


(C12) 


We may straightforwardly define /o = sxQ + h as the posi¬ 
tion of the peak in terms of the figure of merit /, by invoking 
the relation x = s{f — h) which we used to write (Cl). Finally, 
we obtain a set of parameters (/o, A, 7 ), along with a normal¬ 
ization constant yo^ which now all have a direct interpretation 
in terms of features of the modeled distribution (Figure 5c). 
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In summary, they are given in terms of the fitted parameters 
{a 2 ,ai,m) as; 


/o = /r+^^ (-ai+ +8a2m ) ; 


A = (32 


7 tP- 

ZJCq 


- 1/2 


Y= in- 


6x3 > 


(Cl 3a) 
(C13b) 
(C13c) 


recalling that s = ±1 and h in the relation x = s{f — h) are 
fixed by the choice of figure of merit, as given in Table I of 
the main text. The position of the peak is at / = /q. The 
half width of the peak (after it is de-skewed) is given as A (at 
relative height 1 /e), and the factor 7 measures how much the 
peak is skewed towards larger / values (respectively lesser / 
values, if i = — 1 ), with a direct interpretation in terms of the 
horizontal shift of the sides of the peak. 


Appendix D: Confidence regions from the distribntion M(/) 

Here we see how to construct confidence regions from the 
histogram obtained by our method. As explained in the main 
text, regions with high weight in state space may be promoted 
to confidence regions by the method of Christandl and Renner. 

Consider the region with all states p which have at least a 
given value of the figure of merit: 

Rf = {p:f{p)>f} . (Dl) 

The direction of the inequality in (Dl) depends on which 
figures of merit are considered desirable. The direction used 
here reflects the fidelity to a target state, in which case higher 
fidelities are desirable. If, e.g., a proper distance measure such 
as the trace distance is used, the opposite inequality is prefer¬ 
able. 

It is straightforward to see that the weight of the region R f 
in state space according to the measure {p)dp is directly 
given by the weight of the function p (/) over the correspond¬ 
ing range of / values which are included in the region Rf. For 
example, if the figure of merit is the fidelity to a target state, 
the weight a (/) of the region R/ is given by 

«(/)=/ dppB"ip)= f'd/p (/') . (D2) 

JpeRf Jf 

The value of / required for a region Rf to encompass a par¬ 
ticular weight 1 — e/poly(n) is thus given by inverting (D2). 
This may either be done directly from the numerical histogram 
points, or from a fit model. This gives us a region with high 
weight with respect to ps" (p)- 

The method of Christandl and Renner may now be used to 
upgrade these regions to confidence regions. Choose a con¬ 
fidence level 1 — e, and calculate the corresponding poly(n) 
and 5 as given in ref [28]. Recall that a region with weight 
1 — e/poly(n), once enlarged by 5 in purified distance, is a 
confidence region with confidence 1 — e. 


In general, the 5-enlargement can be translated into a cost 
in the corresponding bounding figure of merit /. We’ll de¬ 
rive here this cost for our particular cases of interest of the 
fidelity to a pure reference state, the expectation value of an 
observable and the trace distance to any reference state. Con¬ 
sider first the case where the figure of merit corresponds to the 
trace distance to a reference state PRef, / (p) = D{p, PRef) = 
jIIP — PRefllj. Note the reverse inequality is used in (Dl). 
Then, consider the region Rf+s- Now, we’ll see that Rf+s 
contains the set Rf enlarged by 5 in purified distance. Indeed, 
if p G R/, and a is such that P{p,a) = 2 /l—F^{p ,(t) 5, 

then we may use the triangle inequality, along with the fact 
that D(p,(j) < P(p,(t) [91], to see that D (ff ,pRef) ^ f+5, 
and deduce that a G Rf^s ■ Thus, if is a region with weight 
1 — £/poly (n), then R^ confidence region of confidence 
level at least 1 — e. This construction is depicted graphically 
in Fig. 6 . Of course, the same reasoning applies to the case 
where /(p) = R(p,pRef) is the purified distance to a refer¬ 
ence state. 

Consider also the case where /(p) = tr(plT) is the expec¬ 
tation of an observable W. First, assume that the reverse in¬ 
equality direction is used in (Dl). Then the 5-enlargement 
of Rf is included in the region Rf+„S’ where we’ve assumed 
that the eigenvalues zj of W lie within an interval of size w, 
i.e. w_ ^ Zj ^ w+ and w = w+ — w_, or equivalently, w_lL ^ 
W ^ w+1 and w = w+ — w_. Indeed, assume / (p) ^ / and 
R(p,O’) ^ 5. Then D(p,a) ^ R(p,o) and by properties of 
the trace distance there exists A± ^ 0 such that a — p = A+ — 
A_ and 5 tr(A+ -t-A_) = trA+ = trA_ = D{p,a) ^ 5 [92]. 
Now, /(o) = tr(olT) = f{p) -t-tr(A+VT) — tr(A_VT) ^ f + 
w+ trA+ — trA_ ^ f + w5, and o G Rf+^s- If forward 
direction inequality is used in (Dl) instead of the reverse, then 
the same argument above is easily adapted to show that the 5- 
enlargement of R/- is included in the region Rf_„g, where w 
is defined in the same way. 

Note also that the case of the squared fidelity to a pure ref¬ 
erence state I v/ref) is given by /(p) = (p, = 

(V 6 tef Ip I WRef}’ ^nd is thus the expectation value of the ob¬ 
servable I V 6 tef)(V 6 tef|- More precisely, we have = 1, = 

0 and w = 1 , and the inequality direction used in (D 1 ) encom¬ 
passes larger values for the fidelity in the region; the enlarged 
set to consider is then simply Rf-s- 

Remark that the error bars on the numerical estimate of 
p (/), or on the relevant fit parameters, should a priori be 
propagated to the obtained confidence regions. However since 
p (/) decays exponentially, small errors on the fit paramterers 
will only have a negligible effect on the / required to con¬ 
tain a given weight 1 — e/poly(n) as given by (D2). This is 
just like classical error bars—error bars hardly need their own 
error bars. 

The final confidence regions obtained are still generally un- 
meaningfully large. For example, if we try to construct confi¬ 
dence regions for the two superconducting qubits and choose, 
say, £ = 5%, then we see that £/poly («) ^ 10^37 Yes, that’s 
small. [93] The corresponding / required (see construction in 
Fig. 6 ) is / « 0.85, and we can calculate 5 ~ 0.1. The final 
confidence region comprises then all states with a fidelity to 
the target state in the range [0.75, Ij. This analysis is in itself 
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FIG. 6 . Construction of confidence regions from the histogram of 
the figure of merit. High weight intervals with respect to fl (/) in the 
histogram plots (left plots) correspond to high weight regions in state 
space with respect to (p) (right diagrams). An appropriate en¬ 
largement by a parameter 5 yields confidence regions, a. The case of 
the squared fidelity to a pure target state | V^Ref) figure of merit, i.e. 
/(p) = (p, = (vfiief Ip I Vfiief)• An interval [/, 1] in 

the histogram plot (turquoise region, left) corresponds to a region Rf 
in state space consisting of all states with squared fidelity to | i/tRef) 
higher than / (turquoise region, right). To construct a confidence 
region with a given confidence level 1 — e, first calculate e/poly(n) 
and 8 as given in [28] and choose / such that the turquoise shaded 
area in the histogram plot is 1 — e/poly(«). This means that the cor¬ 
responding region has weight at least 1 — e/poly(n) under p (p). 
Then, the region Rf-S consisting of all states with squared fidelity 
at least / — 5 to | Vfiief) form a confidence region of confidence level 
1 — e (turquoise and orange regions combined) [28]. The true state 
is almost surely in the region Rf^g', equivalently, the true fidelity 
to I V^Ref) is almost surely better than f — 8. Observe that the regions 
constructed this way are linear slices of the quantum state space. This 
is because the squared fidelity to a pure state is linear. Also, the func¬ 
tion p (p) behaves approximately like a Gaussian around the maxi¬ 
mum likelihood estimate Pmle; the illustration of p (p) here corre¬ 
sponds to the case where Pmle coincides with | VtRef)(VtRefl- b. The 
analogous construction applied to the case where the figure of merit 
is the trace distance to the maximum likelihood estimate Pmle- Here 
the regions Rf and R are trace distance balls around Pmle- If f 
has weight at least 1 — e/poly(n), then Rf-^-g is a confidence region 
of confidence level 1 — e. 


not very useful, as it is fair to claim the experiment achieves 
considerably better precision than that (compare with Fig. 3 of 
the main text). The solution we propose is to provide a charac¬ 
terization of the full function p (/) in terms of few parameters, 
from which we know that one can in principle construct con¬ 
fidence regions for any desired confidence level. This is very 
much akin to the eiTor bars reported for a usual classical phys¬ 
ical quantity; these may typically represent one standard devi¬ 
ation of a value which is assumed to be Gaussian distributed. 
A confidence region of high confidence level could then be 
much larger than the reported error bar. For example, a confi¬ 
dence level of one part in a million requires a region size of 5 


standard deviations (or “5 sigma”). 

Appendix E: Application of the method to simulated 
measurements 

We have simulated measurements of individual Pauli oper¬ 
ators on two qubits in the noisy entangled state 

p =0.95 |‘P)('P|-F0.05^ , (El) 

with the pure entangled state |‘P) = (|01) -F(|10)) /V2. Each 
measurement setting consists of a pair of Pauli operators and 
has four outcomes, with a total of 9 measurement settings. 
Each setting was repeated 500 times, resulting in 4500 to¬ 
tal measurement outcomes. Our procedure yields histograms 
corresponding to three different figures of merit (Eig. 7): (a) 
the fidelity to the state I'P), (b) the expectation value of an en¬ 
tanglement witness, and (c) the trace distance to the maximum 
likelihood estimate. The histograms were each generated with 
one random walk instance for each of the 12 CPU cores avail¬ 
able. Each random walk produced 32768 samples, yielding a 
total of 12 X 32768 = 393216 recorded samples. Error bars 
were obtained by binning analysis for each run and combined 
with standard propagation of error bars. This analysis runs 
fast for two qubits and can usually be obtained within minutes 
on usual hardware. 

Eor the figure of merit (a), we have /(p) = 
(p, |‘P)(‘P|) = ('P|p|‘P). This figure of merit is of¬ 
ten used to report the accuracy of experimental preparations 
of quantum states. Here the true value of this figure of merit 
is (p, |‘P)('P|) = 0.9625, given by the “true state” (El) we 
used to simulate measurement outcomes. In (b), the figure 
of merit is /(p) = tr(pW), with the entanglement witness 
W = —1 — Ox C) CTy -F (7y (8> Ox — Oz 0 Oz. The operator W 
is chosen such that tr(pVT) ^ 0 for all states p which are 
not entangled, but also such that tr[|‘P)('P|lT] = 2. Eor the 
last case considered, (c), we first calculate the maximum 
likelihood estimate Pmle, and then define our figure of merit 
as /(p) =Z)(p,Pmle), where D(p,o) := \\\p-a\\i is the 
trace distance. The eigenvalues of the density matrix Pmle 
are (0.000,0.0105,0.0240,0.9655), meaning that the state 
lies on the border of state space. 

Observe that the peak maxima in Eig. 7 do not correspond 
to the values of the maximum likelihood estimate Pmle, even 
though the latter is precisely the point where jib" (p) is max¬ 
imal by definition. This is because of this increasing volume 
factor which shifts the peak. Indeed, at Pmle we can evaluate 
(Pmle, I'/)(’/!) ~ 0.965; the peak maximum in Eigure 7a 
is clearly shifted. 

Let us now apply the fit models to our numerical estimates. 
Eirst, consider the trace distance as figure of merit. The cor¬ 
responding theoretical model is (B8a) with h — 0, as given in 
Table I of the main text; In p (/) « —02/^ — a 1 / + m In / + c. 
The model fits well to the numerical estimate in Eigure 7c 
(red curve). The fit was performed on the logarithm of the 
histogram. We used weights for each point obtained by prop¬ 
agating the error bars as A [Inp] = |(9 (Inp) /(9p| Ap = Ap/p. 
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b. 



j numerical estimate 

— fit approx, model 

— exact Gaussian 


trace distance to MLE estimate 


FIG. 7. Simulated example data analyzed using our procedure, along with the quantum error hars. Two qubits in a noisy entangled state were 
measured with individual Pauli operators, a. Histogram of the squared fidelity to the target state |*P) under the distribution which is relevant 
for the construction of confidence regions using the method of Christandl and Renner [28]. The numerical estimate (blue points) were obtained 
using our procedure based on a Metropolis-Hastings random walk. The theoretical model fits the numerics well (red curve). The quantum 
error bars characterize the distinctive features of the curve, which is interpreted as applying a “skewing” operation on an exact Gaussian (green 
curve). Confidence regions may be constructed from this histogram by choosing regions with states that have a certain minimum fidelity to the 
target state, chosen such that the set has high weight with respect to this distribution. The true value of the squared fidelity of the state from 
which we have simulated measurements is in fact 0.963; the shift is due to the increasing volume factor towards lower fidelity values, b. The 
same analysis is applied to the case of the expectation value of an entanglement witness. The witness is chosen to have positive expectation 
value only for entangled states, with a maximum at the maximally entangled state |'P) where (*? | VT | *P) = 2. c. The same analysis is again 
repeated for the case where the figure of merit is the trace distance to the maximum likelihood estimate. 


Points with obviously huge error bars were excluded from the 
fit. The raw fit for In ji (/) is presented in Fig. 8, along with 
a plot of the residuals. The corresponding fit parameters are 


(with 95% confidence bounds): 

02 = 722.8 

(635.5,810.1) 

fli =319.6 

(305.6,333.6) 

m = 14.09 

(13.82,14.36) 

c = 63.00 

(61.71,64.3). 

The value m is close to M - 

- 1 = 14, which 


(E2) 


would predict from (B6) without the w [/, f2o] term. This indi¬ 
cates that this latter volume factor is approximately constant in 


this case. The quantum error bars can be obtained with (Cl3), 


/o = 0.0377 ; 

A = 0.013; (E3) 

7=0.0014 . 


Our model also fits the numerically obtained histogram in 
Eigures 7a and 7b well. Eor Eigure 7a, the appropriate fit 
model is In^u (/) « —02 (1 —(1 —/) +'«ln (1 —/) + 
c, and for Eigure 7b we have used the model In jj. (/) « 
—(32 (2 —/)^ — ai (2 —/) + min (2 —/)-F c, as specified by 
Table I of the main text. The respective quantum error bars 
are indicated on the corresponding plots in Eig. 7. 
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FIG. 8. Fit of the logarithm of the numerically obtained histogram. 
The log scale allows to better appreciate the quality of the fit. The 
bottom plot shows the residuals of the fit for In fl (/), that is, the dif¬ 
ference between the logarithm of the numerical histogram points and 
the fitted model for the function In/t (/). The low deviation from 
zero for the central points underscores the quality of our model (and 
that the error bars may be overestimated). The points at the sides 
correspond to regions with very low probability and where the nu¬ 
merical estimate is anyway expected to be unreliable. 


Appendix F: Application to experimental data and modeling 
noisy measurements 


corresponding to the real-valued POVM including the rotation 
with Ui are then 


= (Fl) 

(We have ignored here errors in implementing the gate Ui.) 
We could have used these POVM effects directly for each 
measured value for each qubit in the expression for the log- 
likelihood given by Eq. (2) of the main text, however for prac¬ 
tical purposes (to reduce the number of different POVM ef¬ 
fects), we have coarse-grained the values I into 20 different 
bins, yielding the discrete distributions {k) and q[ (k) for 
bin number k. In other words, if the measured value I is in bin 
number k, then the corresponding POVM effect is 

. (F2) 

The joint POVM effect corresponding to combining indi¬ 
vidual measurements on the two qubits is simply given by 
tensoring the two POVM effects. For example, if the value 
I A is measured on qubit A (falling in bin and the value Ib 
is measured on qubit B (which falls in the bin kg), then the 
joint POVM effect is simply 


Q'ij [kAM = Q!i [kA) ® Q!j {ks) , (F3) 

where Ui (resp. Uj) is the rotation applied to qubit A (resp. 
qubit B) before measuring the qubit in the computational ba¬ 
sis. 

We have analyzed the measurement data using the proce¬ 
dure described above. There were in total n = 55677 mea¬ 
surements. The histogram corresponding to the squared fi¬ 
delity to the target Bell state is depicted in Fig. 3 of the main 
text. Our theoretical model (B8b) (with /; = 1) fits the numeri¬ 
cal estimate well. The fit parameters are (with 95% confidence 
bounds). 


As an illustration, we apply our method to analyze mea¬ 
surement data obtained for two superconducting qubits in a 
Bell state prepared using the setup reported in [10]. The data 
were kindly provided by the authors of Ref. [10]. The mea¬ 
surement on an individual qubit is carried out by a transmis¬ 
sion measurement on a resonator coupled to that qubit [78]. 
This measurement yields a random real value I which is dis¬ 
tributed differently whether the qubit is in the |0) state or 
in the |1) state. Single-shot readouts are possible to reason¬ 
able accuracy using a simple threshold, because the two dis¬ 
tributions of I corresponding to |0) and |1) have almost non¬ 
overlapping support [10]. However, we choose to model the 
measurement process more precisely, as our method assumes 
the POVM effects correctly incorporate any noise introduced 
by the measurement device itself. Here, we model the mea¬ 
surement process as a real-valued POVM. A calibration mea¬ 
surement yields the distributions qoil) and q\{I) for trusted 
preparations of the |0) and |1) states respectively. The mea¬ 
surement of the Pauli operator CJ; is performed by applying 
the appropriate unitary gate Ui with high fidelity to bring the 
measurement basis onto the computational basis. The effects 


a2 = 

8511 

(7909,9112) 

ai = 

-476.8 

(-634.8,-318. 

m = 

42.53 

(37.36,47.69) 

c = 

125.4 

(103.5,147.2) . 


From these fit parameters, we finally derive the quantum error 
bars 

/o = 0.934 ; 

A = 0.0086 ; (F5) 

7= 1.4x 10^"^ . 

Appendix G: Comparison with error bars from other methods 

A currently used ad hoc technique for obtaining error bars 
is bootstrapping [22, 24, 27, 41, 79]. In our simulated exper¬ 
iment above, we have performed a simple parametric boot¬ 
strapping analysis for comparison. We have simulated new 
measurement outcomes from Pmle, using the same amount of 
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FIG. 9. Comparison with error bars from existing bootstrapping 
methods. The blue points and red curve reproduce the numerical es¬ 
timate and fit in Figure 7a. We have resampled random measurement 
outcomes from the maximum likelihood estimate and plotted the re¬ 
sulting reconstructed squared fidelities to the target state I'F) (green 
histogram bars). The bias of the curve from our method is due to the 
increasing volume factor in the direction of lower fidelity values. We 
see here that our method yields error bars of the same order of mag¬ 
nitude as bootstrapping. However our error bars are well-justified, 
because they can serve to construct confidence regions for any confi¬ 
dence level. 


measurements and the same settings as for the original simu¬ 
lated experiment. We repeated the procedure many times to 
obtain in total 300 new datasets. For each dataset, we have re¬ 
constructed the corresponding maximum likelihood estimate, 
and determined its squared fidelity to the target state I'P). The 
histogram of these values is presented in Fig. 9, compared 
with the result of our method in Figure 7a. We see that the 
width of the distribution is approximately the same. The bias 
of ~ 1% squared fidelity between the two methods is due to 
the increasing volume factor picked up by our method in the 
direction of decreasing fidelities. Our error bars, however, 
have the robust operational meaning as a means to construct 
confidence regions. 


Appendix H: Discussion of the reliability of the method 

In our work, we provide several levels of reliability state¬ 
ments. First, obviously if /r (/) can be exactly determined, 
then our error analysis is perfectly reliable, assuming the given 
measurement operators are accurate. In practice though, it 
is only possible to approximate jj. (/) with numerical tech¬ 
niques. However these methods are standard and well-tested, 
and come with reliable error estimates [75]. Thus with min¬ 
imal reasonable assumptions an error analysis based on this 
approximation of ji (/) is also reliable. 

Furthermore, we provide an approximate theoretical model 
with only three fit parameters and which explains well the nu¬ 
merical estimate obtained. The quality of the fit over many 


examples studied by the authors not only presents additional 
strong indication that the numerical estimate is faithful, but 
also shows that the result admits a simple representation with 
few parameters. Recall that the numerical method does not 
rely on the assumptions and approximations used to derive the 
theoretical fit model. Also, because of its form the fit model 
is relatively robust to small uncertainties in the fit parameters. 

Our approximate fit model might fail though to describe 
the distribution accurately in some extreme cases, for example 
if too few measurements are taken (however examples with 
n ^ 100 total measurements were well fit). The model is also 
known not to apply to figures of merit such as the fidelity to 
a non-pure state, or more generally, figures of merit which do 
not satisfy (B5). However, the fit in these cases is usually of 
bad quality, especially in the region of the peak. In practice, 
it is sufficient to rely on goodness-of-fit measures or visual 
inspection of the quality of the fit to assert its validity. 


Appendix I: Overview of our software 

We are releasing a software suite which accomplishes our 
procedure in a wide range of settings [83]. The project is com¬ 
posed of a program ready for use, which is built upon a modu¬ 
lar, generic C-n- framework designed for flexibility and speed. 

We expect our program to be directly usable in most ex¬ 
perimental applications. Our program takes as input a list of 
POVM effects £(., which are assumed to be independent, and 
a list of frequencies n^. which indicate how many times each 
corresponding POVM effect was observed. Further inputs in¬ 
clude settings for the histogram range and number of bins, 
which figure of merit to use, parameters of the Metropolis- 
Hastings random walk, the number of times to repeat the ran¬ 
dom walk, the error analysis method, etc. The output of the 
program is the histogram as displayed for example in Fig. 7, 
as well as Fig. 3 of the main text, with corresponding error 
bars. We refer to the project’s hosted location [83] for further 
documentation and detailed information about its usage. The 
histograms presented in this work were all obtained using our 
software. 

This program is itself built upon a generic C-H- framework 
with a collection of tools which may be used to specialize 
our method to more complex setups. We provide for exam¬ 
ple tools to specify the data for a quantum tomography prob¬ 
lem, an implementation of an abstract Metropolis-Hastings 
random walk, an interface to collect statistics during this ran¬ 
dom walk, as well as tools for parallel processing several ran¬ 
dom walk instances. The code is written using a technique 
called Ch-h- template metaprogramming [94], which allows to 
write generic code which is flexible and reusable, but which at 
compile-time is translated into highly optimized low-level ma¬ 
chine instructions. Our project relies on the Eigen and Boost 
libraries [95, 96], in particular for linear algebra calculations. 
These libraries also make extensive use of this technique. 

Some tasks are not covered by our program. If required 
by the figure of merit, finding the maximum likelihood esti¬ 
mate Pmle can be accomplished by minimizing the loglike- 
lihood A (p). Since this function is convex, the solution can 
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be found efficiently. In most of our examples, we used CVX, 
a MATLAB package for specifying and solving convex prob¬ 
lems [97, 98]. Also, in order to determine the fit parameters 
corresponding to the histogram, we resorted to MATLAB’s 
curve fitting toolbox. A MATLAB script is provided along 
with the software to ease this task, and to calculate the quan¬ 
tum error bars. 


Our program is currently limited to POVM effects with a 
product structure. However, this is generically the case, even 
e.g. for adaptive tomography [99-101]. We have successfully 
used our code to analyze simulated measurements of Pauli op¬ 
erators of up to at least 5 qubits on our hardware, and we ex¬ 
pect further improvements will increase this limit. 
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